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Classical fields approximation to cold weakly interacting bosons allows for a unified treatment of 
condensed and uncondensed parts of the system. Until now, however, the quantitative predictions 
were limited by a dependence of the results on a grid chosen for numerical implementation of the 
method. In this paper we propose replacing this unphysical ambiguity by an additional postulate: 
the temperature of the gas at thermal equilibrium should be the same as that of an ideal Bose gas 
with the same fraction of condensed atoms. As it turns-out, with this additional assumption, nearly 
all atoms are within the classical fields, thus the method applies to the whole system. 

PACS numbers: 03.75.Hh 
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I. INTRODUCTION 

The experimental realization of a Bose-Einstein con- 
densate in dilute alkali atomic gases [E H; Q has opened 
a novel possibilities to explore the behavior of large quan- 
tum systems. A study of an interplay between classical 
and quantum regimes is now available, since the experi- 
ments are carried out at very low, yet non-zero tempera- 
tures. This is an important issue, because many of these 
results depend on the temperature. 

The atomic BEG is a many-body quantum interacting 
system, thus it is not easy to build a usable dynami- 
cal theory, that takes thermal effects into consideration. 
One group of models addresses this task in a two-gas ap- 
proach, by explicitly dividingthe system into the conden- 
sate and the thermal cloud 4, 5, 'G\i\- Both subsystems 
evolve with different equations and they influence each 
other by a mean-field interaction. Such an arbitrary sep- 
aration does not however allow for a deeper insight into 
the dynamics of the BEG. 

The second approach, explored in this paper, is the 
classical fields approximation (GFA, or its canonical ver- 
sion known as the Wigncr function method) 8, 9, 10, llj. 
It is a mean-field approximation of the quantum field 
theory of interacting bosons, in which one identifies the 
"classical modes" and evolves them using the Gross- 
Pitaevskii equation (GPE) on a grid. The need to choose 
the finite number of modes for simulation introduces a 
cut-off parameter into this method. The approximation 
does not allow to describe quantum correlations, but it 
lets one observe the emergence of the BEG from the pure 
Hamiltonian dynamics. This is remarkable as condensed 
and non-condensed atoms are described in the same way 
in the classical fields approximation, and the interaction 
and the observation process allows to distinguish between 
the two phases. 



In this paper we present a version of the classical fields 
approximation in a 3D box potential without free param- 
eters. This goal is achieved by introducing a transforma- 
tion from numerical control parameters: the energy per 
particle E, an interaction strength g and a grid size Ugrid, 
to physical control parameters: the number of particles 
N, the temperature T and the scattering length a. Such 
transformation is not uniquely defined; in fact, it consti- 
tutes the major difficulty in interpretation of the classical 
fields approximation. The main problem in applying the 
GFA has been to properly tune the grid size ngrid to 
the temperature and the number of particles. In previ- 
ous applications authors used different approaches. One 
of them is to fix the population of the highest momen- 
tum mode to some arbitrary value and calculate number 
of particles a posteriori [13, [H, U^ uM ■ The other one 
[3, [l3 is t° set the numerical grid and the number of 
atoms on the grid, while treating the classical fields as 
just a small subset of a bigger system. 

In this paper we suggest yet another method of assign- 
ing the physical parameters to a given numerical simu- 
lation. The main idea is to set the temperature of the 
system to the value of an ideal gas with the same conden- 
sate fraction. This is of course an approximation, since 
it is known that the critical temperature depends on the 
interaction strength. However, in the weakly interacting 
regime a shift of the critical temperature is very small. 

The paper is organized as follows. In Section ^] we 
introduce the classical fields approximation. We identify 
numerical control parameters that determine the prop- 
erties of the equilibrium of the system. Next we review 
eigenmodes of the system - the generalized Bogoliubov 
quasiparticles. 

In Section IIIII we present an analysis of the thermal 
equilibrium of the system. We observe a phase transition 
as an abrupt change of properties of the system and show 



that the equipartition of energies in eigenmodes of the 
system occurs in a condensed phase only. We determine 
the dimensionless teniperature per particle r, which is 
proportional to T /N [l2. Il5| . We also explain a pitfall 
that one encounters when trying to define the method 
via setting unphysical parameters: the grid size or the 
population cut-off. 

In Section IIVI we compare the scaling properties of 
the temperature of a system obtained from equipartition 
of energies with the temperature of the ideal Bose gas. 
From this we derive a formula for the number of particles 
and obtain a set of physical control parameters (iV, T, 
a, L) without referring to the grid size. Next, we per- 
form an approximate check of consistency by comparing 
the numerical kinetic energy per particle with a result 
for an ideal gas and find an agreement within 40%. Fi- 
nally, we confront the calculated population cut-off with 
assumptions of the classical fields approximation. 

In Section we present some conclusions. 



II. THE METHOD. NORMAL MODES AND 
THEIR ENERGIES 



Let us consider a gas of N identical bosons trapped 
within a 3D box potential of length L with periodic 
boundary conditions. Wc assume that atoms interact via 
a contact potential V{r — r') = ATrh'^aS{r — r')/m, where 
a is the s-wave scattering length, known to be adequate 
at low temperatures. 

The second-quantized Hamiltonian reads: 
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I £x (^t|^^) + ^^^ / d'x (^t^t^^) 
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where ^ is a bosonic field operator satisfying equal time 
bosonic commutation relation ['^{r,t),'i'^{r',t)] = S{r — 
r'). The Hcisenberg equation for ^ resulting from this 
Hamiltonian is of the form: 

^t,j, ^ (2.2) 



ihdt'H = - 
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The symmetry of the box with periodic boundary con- 
ditions sets a natural basis of plane waves with quantized 
momentum p = 2'Kh'k/L, thus it is convenient to expand 
the field operator: 
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(2.3) 



Annihilation operators at destroy particle in mode (fc) 
and satisfy a commutation relation [at, aj ^,] = 1. 

Using the above decomposition in Eq. i|2.2|l we get a 
set of nonlinear operator equations for annihilation oper- 
ators of plane wave modes: 



„ . / s 2Tr'^h 2 * / N 47rfta ■r-^ 

Otai^[t) = -I — —^ k ai^[t)~'i — ^v > a' 



mL^ 



qi,q2 



qjOq^ak+qj-q^ 

(2.4) 



A full solution of these equations is not known. How- 
ever we can simplify this problem by applying an approxi- 
mation which is an extension of the Bogoliubov approach. 
We extract from the field operator all modes that are oc- 
cupied by sufficiently large number of atoms, such that 
one can justify neglecting their quantum nature, and re- 
place their annihilation operators with c-numbers: 



flk 
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(2.5) 



This approximation has been extensively used in de- 
scription of a multimode laser light in quantum optics; 
here it corresponds to an assumption that these " classical 
modes" are coherent. The square of the modulus of the 
amplitude |akP gives the fraction of atoms that occupy 
the mode k. 

Our set of classical modes corresponds to momenta 
forming a 3D cubic grid with Ug points in each direc- 
tion. The size of this lattice is defined by the "classical 
mode" of the largest momentum. This way we introduce 
maximal momentum into the problem. This momentum 
cut-off is the only additional parameter in the method. 

Substituting operators for " classical" modes with their 
corresponding classical amplitudes in equation 12.41 and 
neglecting all remaining operators we get: 



, , 2Trh , , , , AnhaN ^^ 

Otauit) = -« — 77 fc ak(i)-« 



ki, k2 

(2.6) 

We introduce for convenience a unit of energy e = 

Att'^H'^ /{mL'^) and a corresponding unit of time: h/e. 

The Eq. 12. 61 can be rewritten in these units in 

terms of a dimensionless mean-field wave function ^ = 



(2.1) Ek^ke-^-'^'-/^ 



as: 



idtip = -yV + .glV'I^V', 



(2.7) 



where g — aN/nL is the interaction strength and the 
unit for ip is L~^^^. It is the celebrated Gross-Pitaevskii 
equation on a grid and can be effectively solved by means 
of split operator method using Fast Fourier Transform 
(EFT). 

An important feature of Eg 12. 71 is that for almost all 
initial states ■tp{t = 0) for a given grid, scattering length 
a and a box size L the evolution leads, after a tran- 
sient period of time, to a thermodynamically steady state 
that depends only on the energy per particle [ll|, [lj|. 
This thermalization is illustrated in Fig. ^ where we 
present the time-evolution of the population of the zero- 
momentum mode |ao.o,o(i)P- Thus the number of con- 
trol parameters for equilibrium states is reduced to four: 
ngrid, E, L, g. 

The wave function ip stands for both the condensate 
and the thermal cloud, as opposed to standard inter- 
pretation where ip stands for a pure BEC and only the 
ground state of the GP equation is considered. In fact 
we have not made any distinction between condensed and 
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FIG. 1: Time-evoluiton of a zero momentum component 
|Qo,o,o(i)|'^- A thermalization is clearly visible, as well as rapid 
fluctuations due to non-linear dynamics. The data have been 
obtained on a grid 36x36x36 and g — f8.5f5. 



thermal atoms. Instead we use a classic criterion of On- 
sager and Penrose 16] to identify condensate fraction as 
a dominant eigenvalue of the single particle density ma- 
trix. 

A careful reader will notice that p{r,r') — ip{r)*ip{r') 
is a pure state and has only one eigenvalue. However, 
typical measurements of BEC involve optical techniques, 
with the exposure time. At, varying from a few up to 
hundreds of milliseconds. On the other hand, the evo- 
lution of ip is very rapid and irregular (note fluctuations 
in Fig. ^1 due to a very short time scales of a nonlinear 
many-body dynamics. An observation leads to a coarse 
graining and the measured density matrix is of the form: 



Paver yi^) — 



At 



t+At/2 



t-At/2 



dr ^*{r',T)^P{r,T) (2.8) 



This time averaging procedure destroys the purity of 
the state by destroying coherence between eigenmodes. 
It also smoothes out the irregularities so that time- 
averaged density profiles greatly resemble experimental 
photographs [l4| . This way we can interpret our iso- 
lated system as a mixed state - a single realization of a 
quantum system with the properties of a measurement 
process taken into account, which is a case similar to 
experiments. 

In general, to obtain eigenmodes of the system one 
needs to diagonalize the time-averaged density matrix 
Paver = Sk "■k0k(^')^k(?') , whcrc (j)k{r) are eigenvecors 
of Paver {r,r') (scc ll] for application of CFA to a har- 
monic potential). However, the diagonalization is not 
required for a discussed 3D box potential with periodic 
boundary conditions due to its symmetry [I2j |: nk are 
just time averaged squares of amplitudes of mode k, 
rik =< |ak(OP >; ^ii'i physically they are equal to rela- 
tive populations of eigenmodes. 

This way both dynamic and thermodynamic properties 
can be simultaneously studied |12.ll3lll^ . Still the main 
problem of the CFA and the goal of this work is to assign 
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FIG. 2: Typical frequency spectra of plane-waves modes. 
The condensate mode (bottom-left) has a clearly defined en- 
ergy p, while spectra of excited modes are broadened due to 
interactions. Note that one group of peaks is suppressed in 
(0, 0, 5) mode - the Bogoliubov transformation is practically a 
unity in this case. The parameters of the plot are: rigrid = 32, 
g = 18.515, no = 86.2%. 



results to meaningful physical parameters: the number 
of particles N and the temperature T. To achieve this we 
now review normal modes of the system - the generalized 
Bogoliubov quasiparticles, which have been thoroughly 
studied in ilj| . We do this for the sake of completeness 
of this paper, as we need it to present equipartition of en- 
ergies in Section lHll which in turn allows us to determine 
T and N from numerical control parameters. 

Typical frequency spectra of time-evolved amplitudes 
Ok are depicted in Fig. El For excited modes, they 
consist of two groups of peaks centered at values /i ± ek, 
where ek is given by the gapless Bogoliubov-like formula: 



Ck 



fc2 
(Tr+5"o)^ - (5"-o)^ 



(2.9) 



and p is the energy of the condensate mode. For large 
momentum this formula reduces simply to fc^/2. With in- 
creasing momentum the left group of peaks is suppressed, 
such that for high momentum modes only one group of 
peaks remains. 

The approximate equations for excited modes couple 
amplitudes ak with a^^. (TJ. Provided that spectra of 
these modes consist of two peaks instead of two groups 
of peaks, we can perform Bogoliubov transformation to 
obtain a quasiparticle amplitude (5k oscillating with only 
one frequency p+e^. In our case we deal with two groups 
of peaks, thus we need to generalize this reasoning. We 
decompose ak and a_k into corresponding pairs of spikes 
with the frequencies {p+u) and (p — to). By applying the 
Bogoliubov transformation to each pair of peaks we arrive 
at the (p+Lu) component of the frequency spectrum of the 
quasiparticle amplitude: Su_{p + u). Their squared am- 



plitudes are occupations of quasiparticles n^" 
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FIG. 3: Frequency spectra of quasiparticles - the eigenstates 
of the system. They have only one central energy, similar to 
high-momentum plane-wave modes (compare with bottom- 
right picture in Fig. 12.91 . The broadening of the spectrum 
determines the lifetime of a quasiparticle. The data was ob- 
tained on a grid Ugrid ~ 32 and for the interaction strength 
g — 18.515. The condensate population is 86.2%. 



These quasiparticles represent the eigenmodes of the sys- 
tem with their energies being the central frequency of 
their spectra (Fig. |3), equal to /x -f et- The detailed for- 
mulas for this procedure are presented in the Appendix. 
At the end of this section, let us introduce the cut-off 
population Numax to conveniently confront our method 
with literature and the assumption of the classical fields 
approximation. We define it as a mean population of 
" edge" classical modes with momenta greater or equal to 
nhngrid/L. Such modes are included in our calculations 
because we work with rectangular grid. The population 
Numax is the population cut-off that separates classical 
(mean-field) and quantum (neglected) modes and quanti- 
tatively specifies the term " sufficiently large population" . 
It can be used as an alternative numerical control param- 
eter instead of the grid size Ugrid [H, U^ , as Numax is 
defined by the grid for given L, g and E. 



III. A TRANSITION TO BEC. A 
DISTRIBUTION OF ENERGIES. 

In this section we analyze properties of the system in 
thermal equilibrium above and below the critical energy. 
We use the equipartition of energies in eigenmodes of the 
system to determine the dimensionless temperature per 
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FIG. 4: A probability distribution of occupations of excited 
modes. Plots in insets are in logarithmic scale to outline the 
observed exponential trend. The parameters are: Ugrid ~ 32, 
no = 86.2% and g = 18.515. 



particle r. Finally, we discuss in detail the inconsistency 
emerging from the free choice of the grid for numerical 
calculations. 

First, we inspect the behavior of excited modes in the 
state of equilibrium. We define pw{n) to be a probability 
of finding the population of quasiparticle n^"""^ equal to 
n. It is calculated by counting all events where |(5kP is 
close to n during a simulation. Results for pk are pre- 
sented in Fig. ^ They reveal roughly exponential dis- 
tribution of populations of thermal modes and provide a 
strong argument that we observe the thermalization in 
the system. Thus one can treat a single excited mode as 
being in thermal equilibrium with a reservoir consisting 
of the rest of the system. 

For comparison, similar calculations of po,o.o('^) per- 
formed for the condensate mode ao,o,o reveal a phase 
transition into BEC ([H IH 13). Below the critical 
energy the distribution is peaked around a non-zero con- 
densate population (Fig. ISji), while above E^ it remains 
exponential (Fig.[Sl3). 

Let us now consider a distribution of populations of 
thermal modes n^. We have found that this feature dif- 
fers for condensed and uncondensed systems (see Fig. ^ 
and defines the critical energy more accurately than the 
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FIG. 5: A probability distribution of occupations of the con- 
densate mode: (a) below Ec {g = 18.515, ng=48, no — 7.2%), 
(b) above Ec (rig = 50, g = 18.515, E = 300). 



change in po,o,o("-)- Below the Ec the distribution is pro- 
portional to k~^ and this agrees with the low tempera- 
ture limit for the population distribution of an ideal Bose 
gas. Above the critical energy the distribution abruptly 
changes to exp(— fc). 

The distribution of populations is closely related to the 
equipartition of energies occurring in eigenmodes of the 
system: 



Nrikek = ksT. 



(3.1) 



This important concept, introduced to the CFA in [13 
and expanded in |l2l |. allows one to easily determine the 
dimensionless temperature per particle t - the essential 
parameter for our method. 

We find that the equipartition occurs only below the 
critical energy, thus limiting its application to calculate 
the temperature to below the critical temperature Tc- 
One can notice this easily for high-momentum modes 
(see Fig. [7^), where e^ = h^k'^/2m and the equipartition 
reduces to n^fc^ = const. Clearly, the exponential popu- 
lation distribution occurring above the critical energy Ec 
cannot satisfy this relation (Fig. [7]^) . 

Before we make use of the equipartition to obtain r, 
let us present yet another argument supporting the qual- 
ity of observed equilibrium - the fluctuations of energy, 
visible in Fig. [7Ji. A system which is truly thermalized 
experiences such fluctuations, but they vanish with the 
averaging time as Ai~^/^, according to the Central Limit 
Theorem. Indeed, the agreement of results with the the- 
orem is excellent, as can be seen in Fig.|Sl 

Now we can determine the dimensionless temperature 
per particle r = ksT /eN . The equipartition can be writ- 
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FIG. 6: Populations of thermal modes versus their momen- 
tum: below (a) and above the critical energy (b) . Parameters 
are: (a) Ugrid = 54, .B = 215, no = 0.3%, g = 18.515, (b) 
i^grid = 50, E = 225, g — 18.515. Note how a very small 
fraction of condensed atoms (0.3%) can influence the whole 
system. 



ten in a form independent of A'^ and L as: 
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71k 



(3.2) 



For different condensate fractions uq this relation is de- 
picted in Figure |^ The parameter r is a slope of the 
best linear fit to numerical data [Tjl ■ It depends on three 
numerical parameters: the energy per particle E, the in- 
teraction g and the grid size Ugrid, such that if g is set, 
then a pair (no, t) is unequivocally determined by E and 

^grid- 

There are two common ways of specifying the classi- 
cal fields approximation in the literature. One of them 
[12. Il4| defines quantitatively the assumption of "suf- 
ficiently" large population of modes treated classically 
by arbitrarily setting the value of the population cut-off 
NiT'max- The other method 'y, '15] assumes that numer- 
ical results represent only a fraction of a bigger system; 
the population cut-off is large (i.e. greater than 15) and 
external atoms are approximated by an ideal gas. 

The major concern of the CFA is the dependence of 
results on the particular choice of the grid, as it does not 
have any physical meaning. It is illustrated in Fig. 1101 
which depicts the condensate fraction no versus the tem- 
perature T = teN/L obtained for various grids Ugrid and 
energies per particle E, while g and L are set. We have 
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FIG. 7: The energy confined in excited modes (ekWk) versus 
the momentum for g — 18.515: (a) below Ec, Ug^id ~ 48, 
E = 163.6 and (b) above Ec, Ugrid = 50, £ = 300. Only 
below the critical temperature (no = 7.2%) one observes the 
equipartition of energy. Note that the condensate mode does 
not satisfy the equipartition. 
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FIG. 8: Fluctuations of an inverse dimensionless temperature 
per particle (3 = r~^ as a function of the number of averaging 
time-steps. The data is obtained on a grid Ugrid ~ 48 and the 
condensate population is no — 7.2%. 



FIG. 9: The energies of quasiparticles versus their inverse 
populations depicted for a condensed system (no = 7.2%, 
ngrid = 48, E = 163.6, g = 18.515). A slope of the linear fit 
to the data is the dimensionless temperature per particle r. 



One can clearly see that by changing the numerical 
grid it is possible to obtain a relatively broad range of 
temperatures T (and t) for a single condensate fraction 
no. The grid that reproduces the result for the ideal Bose 
gas at one value of no, fails to do so for different conden- 
sate fractions. On the other hand there exists an optimal 
value of the population cut-off Nrimax (ca. 0.6 — 0.7) for 
which the results match the ideal Bose gas both for low 
temperatures or close to T^. If Nrimax differs from the 
optimal value not too much, the results follow a curve 
which is shifted with respect to the ideal gas (see results 
in [12| , where the population cut-off is arbitrarily set to 
1). Such results are still reliable, although their accuracy 
is reduced by several tens of per cents. However, the 
optimal value of the population cut-off Nrimax changes 
with the interaction strength g, so that there is no uni- 
versal value that satisfies all control parameters. And as 
its physical interpretation can be qualitative at most, the 
population cut-off is still a free parameter of the method, 
supplementary to the grid size Ugrid- 

To summarize, we have analyzed the properties of the 
thermal equilibrium of the system and shown difficulties 
arising from the existence of artificial parameter (Nrimax 
or Ugrid) in the classical fields approximation. We have 
based the quality of results on the comparison of the 
temperature calculated from equipartition with the cor- 
responding temperature of the ideal Bose gas. In the 
next section we build on this concept to obtain a method 
which is free of the presented vulnerabilities and yet de- 
scribes the whole gas within the classical fields in a simple 
and consistent way. 



assumed N to be constant (235000) so that we can deter- 
mine T and the population cut-off Numax- A solid curve 
represents a relation no{T) of the ideal Bose gas with 
the same N and L. For interacting Bose gas this relation 
would be only slightly modified, as for small scattering 
length corrections to the critical temperature are known 
to be small Ullllli^. 



IV. ELIMINATING THE CUT-OFF 

We base this approach on an assumption that numer- 
ical results (no, r) represent a certain entire Bose gas 
with the condensate fraction no. We impose a condi- 
tion that the temperature obtained from the equiparti- 
tion T = '^\^1-i^ T matches the temperature of the corre- 
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FIG. 10: The condensate fraction versus the temperature 
for various numerical grids. The interaction strength g is 
18.515 and the number of atoms A*' is set to 235000. Num- 
bers describing the points in the figure are the population 
cut-off Nrimax (left numbers) and the grid size Ugrid (numbers 
in brackets). The population cut-off is known only approxi- 
mately due to large fluctuations in highly excited modes. The 
optimal value of the cut-off, for which the numerical results 
agree with the corresponding ideal Bose gas is ca. 0.6 — 0.7 
in this case. 



sponding Bose gas, approximated with a formula for an 
ideal Bose gas: 
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(no,iV,L) = l^r'/3(3/2)(l-no)2/3 



mk 



7V2/3 
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(4.1) 
where C is the Riemann function. 

To justify this approximation, note that for low scat- 
tering lengths the critical temperature of an interacting 
Bose gas is very close to that of of an ideal Bose gas, 
although there is a significant disagreement between var- 
ious authors over its precise value (see for instance [13 
for calculations of the corrections). For typical systems 
with the atomic density of 1.8 IO^'^cttt,^'^ and the scat- 
tering length a = 5.8nm the correction to Tc is about 
4.5%. 

Because formulas for both temperatures scale differ- 
ently with the number of particles, we can always satisfy 
our assumption by taking: 
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(1 - no? 

(2(1) ^3^3 



(4.2) 



This way we determine the number of particles N and 
the temperature T of the system from the condition that 
has a clearly defined physical meaning. The remaining 
physical control parameter - the scattering length a can 
be now easily calculated from the value of the interaction 
strength g = aN /-kL used in numerical calculations. 

The main advantage of this approach is that we do not 
set arbitrarily any artificial parameter, because N is im- 
plicitly dependent on the grid size Ugrid- Thus problems 
arising from a freedom of choice of the numerical grid, 
such as outlined in Section [llll do not occur anymore. 



We perform a rough check of consistency of our result 
by comparing the kinetic energy per particle obtained 
from the simulation with the kinetic energy per particle 
of an ideal Bose gas, which is given by the formula: 



E, 



kin 



N 



3 3 i, 5 , ,^a 



r 



(4.3) 



The results show an agreement within 40% , what is 
acceptable due to the very approximate nature of this 
comparison. 

Next, we calculate the cut-off population Nn^rnax- It 
represents the optimal value of the cut-off for chosen 
condensate fraction tiq and interaction strength g. ft 
is remarkable that, even though we change the grid and 
adjust the energy per particle such that no remains con- 
stant, our procedure alters the number of particles in such 
a way, that the cut-off remains relatively fixed (approx. 
0.7 for g — 18.515). In general, it depends only on the 
interaction strength, and this agrees with our previous 
observations (see Chapter IIII|) . 

The relative invariance of the population cut-off with 
respect to the grid size rigrid and the condensate fraction 
no can be understood by considering a simplified model 
of a system with equipartition of energies in plane waves 
instead of proper eigenmodes. Such system (considered 
earlier in Q) corresponds to the case of g = 0, thus it 
cannot be obtained within the CFA, which relies on inter- 
actions for thcrmalization. We calculate the number of 
particles again by comparing the dimensionless tempera- 
ture, given by the relation r = (1 — no)/(X]k w) ^^*^ *^^ 
temperature of an ideal Bose gas. The cut-off population 
for this system is given by the formula: 



Nrir, 



(^k 



j-12 y. 



U\>Tr7igr-id/L TP 



^'CHD E|ki>. 



r-id/L 



1 



(4.4) 



Similar to the interacting case, this population cut-off 
does not depend on the condensate fraction no. For the 
grid size Ugrid — 16 it equals 0.69, while at ngrid ~ 128 
it is 0.80. This change in the population cut-off can be 
attributed to the box shape of the grid. Indeed, if we 
assume a spherical grid and approximate summations 
in Eq. 14.41 with integrals, then the result {Nrimax = 
IGtt^^C 2(3/2) = 0.75) is independent also from the grid 
size. 

Note that the presented approach does not depend on 
the value of the population cut-off (or any other artificial 
parameter). This parameter serves only for comparison 
with the assumption of " sufficiently large" populations of 
classical modes. It can be surprising that the optimal cut- 
off population comes out so small, but on the other hand 
it shows that almost whole gas is confined in the classical 
modes. An estimate based upon approximating external 
modes with an ideal gas and gluing both gases with a 
single temperature and the same cut-off population on 
the border yields about 97 — 99% atoms inside the CFA 
grid. This result is consistent with our assumption that 
we describe the entire system. 
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At the end, we also discuss two drawbacks of the pre- 
sented approach. First, it does not allow to study the 
shift of critical temperature upon the scattering length, 
as we fix the temperature to the value of an ideal gas. 
The second limitation is that one can apply this method 
only below the Tc^ as we use equipartition of energy which 
occurs in a condensed phase only. However, we are very 
uncertain if the CFA works above Tc (at least in 3D box), 
as the observed exponential population distribution of 
excited modes does not agree with the distribution of an 
ideal Bose gas. Also, the CFA has not been constructed 
for this purpose. 



V. CONCLUSIONS 

We have presented the version of the Classical Fields 
Approximation without free parameters, which finally al- 
lows one to easily and consistently present results versus 
real- world parameters, like the temperature and the num- 
ber of particles. Most important, the accuracy of these 
results is unbiased by the choice of the numerical grid - a 
major concern of former applications of the CFA method, 
as we rely on the condition with a clear physical meaning 
rather than arbitrarily set some artificial parameter. 

We have calculated the optimal population cut-off lim- 
iting the classical and quantum regimes in the CFA and 
found that it almost does not depend on the condensate 
fraction or the grid size. We explain this behavior with 
a simplified model of a non-interacting gas on a grid and 
find considerable agreement. As a conclusion we validate 
the previous use of the fixed cut-off classical fields approx- 
imation [121, although we point to the hardly controlled 
accuracy of such results. 

Moreover we have shown that one can describe the en- 
tire gas within the classical fields, as opposed to an ap- 
proach presented in [3, where the CFA is considered to 
describe only a small part of the whole system. However, 
from the obtained value of the population cut-off we con- 
clude that the classical fields approximation is exploited 
to its limits of applicability. 

We have also presented the detailed analysis of the 
thermal equilibrium of the system. We have reported 
on the rapid change in the distribution of populations 
of the excited modes which marks the phase transition 
into BEC in the classical fields approximation. We have 
shown that the equipartition of energies in eigenmodes 
of the system occurs in a condensed phase only, which 
limits the use of our method to temperatures below the 
critical one. 
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APPENDIX: A BOGOLUBOV 
TRANSFORMATION 

We present here the derivation of formulas for the gen- 
eralized Bogoliubov quasiparticles, which have been de- 
scribed in Section^ The observed spectra of amplitudes 
of excited modes consist of two groups of peaks (see Fig. 
|2Jl, thus we describe the amplitude of mode k as: 



ttk 



[t) = V /3k(c.)e-'(^-")* + 7k(cc')e"'(^+")*, (A.l) 



where the coefficients (3\s_{uj) and 7k (w) are obtained from 
simulation. 

As the approximate equations for excited modes cou- 
ple amplitudes ctk with al^. |12J| . we can identify the 
corresponding pairs of spikes Q;k(M + ^) ^iid a-k(/^ — ^)'- 
they are 7k(w) and /3_k(w), respectively. We apply a Bo- 
goliubov transform to each such pair to obtain a single 
spike having only one frequency. In order to do so, we 
introduce a new quasiparticle amplitude (5k, defined as 
follows: 

6v,{^i + Lo) = Ui,{uj)a±i,{fi + uj) + e-2vVk(cj)alk(/i - uj), 

(A.2) 
where Uk^cu) and Vk(a;) satisfy the condition: 

|C/±k(^)P-|l^±k(^)P = l (A.3) 

Substituting lA. l| into IA.21 we get: 

(t/k(c^) 7k(^) + ^k(c^) /31k(^))e"'^^+"^*- 

(A.4) 

We want (5k (/i + uj) to be a component of a quasipar- 
ticle with only single positive energy, thus we impose a 
condition: 

Ukiiu) /3k(c^) + Vk(c^) 7lk(^) = 0. (A.5) 

From eq. I A. 31 and I A. 51 we obtain: 

l7-k(c^)| 



Uuiuj) 



Vk{L^) 



l/?k(u;)| 
Vl7-k(c.)P-|/3k(^)P 



-iArgiMio)) 



oJ-4rg(7_k(ij)) 



(A.6) 



(A.7) 



Finally we arrive at the fi + uj component of the quasi- 
particle amplitude composed of modes k and — k: 



Sk{^i + uj) 



-i{fi-\-Lo)t 



Vl7-k(u;)P~|/3k(c.)P 



(l7k('^)7-k(w)|e 
l/3k(w)/3-k(t^)|e'(^''»''^-''<'^^^-^'^»'''-''^'^))) 



Spectra of these excitations are centered on a single 
value (Fig. ^, similar to high momentuni plane- wave 
modes. They can be regarded as normal modes of the 
system oscillating with approximately single frequency, 



while the width of the spectrum is related to the life-time 
of a quasiparticle [i3 • 
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